TK!2
=v
=u
=l
#1
:n
a0
:v
'a01
'a02
'a03
'a04
'a05
'a06
'a07
'a08
'a09
'a010
E
#2
:n
a01
:v
74
-25
66
-11
-96
22
-40
-2
-22
-3
E
#3
:n
a02
:v
58
54
-50
36
-21
-34
23
-64
44
-58
E
#4
:n
a03
:v
-37
-100
-41
-59
-99
52
75
42
99
87
E
#5
:n
a04
:v
81
-68
79
87
-39
55
-62
-45
-2
18
E
#6
:n
a05
:v
-86
17
-58
-71
26
-38
35
-53
90
33
E
#7
:n
a06
:v
-25
81
-93
-78
-11
48
57
80
14
-90
E
#8
:n
a07
:v
-56
95
5
43
68
40
43
-34
-89
20
E
#9
:n
a08
:v
86
46
22
8
25
42
-87
27
-35
99
E
#10
:n
a09
:v
-44
-79
-79
-94
83
68
58
-76
-67
-84
E
#11
:n
a010
:v
-84
-24
56
25
63
22
-55
-24
47
-26
E
#12
:n
b0
:v
-694
-684
1843
33
353
44
329
935
-1637
169
E
#13
:n
a
:v
'a1
'a2
'a3
'a4
'a5
'a6
'a7
'a8
'a9
'a10
E
#14
:n
x
:v
E
=f
#1
:n
Mcopy
:c
Copy matrix (i.e. subordinate lists)
:t
Procedure
:i
M1,M2
#2
:n
Mdelete
:c
Delete matrix (i.e. subordinate lists)
:t
Procedure
:i
M
#3
:n
gelpiv1
:c
Gauss Elimination with Scaled Partial Pivoting
:t
Procedure
:i
n,a0,b0,a,b,x
#4
:n
gelpiv2
:c
Gauss Elimination with SPP - iterative improvement
:t
Procedure
:i
n,a0,b0,a,b,x
#1
:n
Mcopy
:z
" The procedure copies subordinate lists of the masterlist M1 into 
" subordinate lists of the master list M2.

for i=1 to length(M1)
   call listcopy(M1[i],M2[i])   " error indicates no listname in ith
                                " element of either M1 or M2
next i
/E
E
#2
:n
Mdelete
:z
" The procedure deletes subordinate lists of the master list  M  from the
" List Sheet. The master list remains intact.

for i=1 to length(M)
   call delete(M[i])    " error indicates no listname in ith element of M
next i
/E
E
#3
:n
gelpiv1
:z
" Notation:  n     - number of equations
"            a0,b0 - lists with given matrix and RHS coefficients
"            a,b   - working lists with matrix and RHS coefficients
"            x     - list with solution result
"            scale - list with scaling coefficients
"            idx   - list with row pointers
" Related functions:  Mcopy, Mdelete

" Description: The procedure 'gelpiv1' uses one of the most popular and
" robust algoritms for solving systems of linear equations.  The scaling and
" pivoting controls the effect of round-off errors.  The singularity of
" a matrix is tested against a 'floating point zero' which takes into account
" the effect of the number of digits, and of the size of the problem.

" Use: The matrix coefficients have to be stored in lists row by row. 
" The names of the given and working lists must be specified as input
" arguments. Sufficiently long lists of row listnames for the working lists
" must be set up prior to calling procedure gelpiv1.



call Mcopy(a0,a)          " copying matrix and RHS coefficients
call listcopy(b0,b)       " from input lists to worklists

zero:= n*1e-16            " floating point zero threshold

for i=1 to n              " generating scale coefficients for pivoting
   'idx[i]:= i
   big:= max(max(a[i]),-min(a[i]))
   if big<zero then error:= matrix_with_zero_row
   'scale[i]:= 1/big
next i

for k=1 to n-1              " Gauss elimination
   big:= 0    
   for i=k to n             " selecting pivot
   idxi:= 'idx[i]
   new:= abs(a[idxi][k])*'scale[idxi]
   if new<=big goto A
     big:=new
     idxp:= i
A: next i
   if big<zero then error:= singular_matrix
   if idxp>k then ('idx[k],'idx[idxp]):= ('idx[idxp],'idx[k])
   idxk:= 'idx[k]
   pivot:= a[idxk][k]
   for i=k+1 to n           " matrix transformation
      idxi:= 'idx[i]
      aop:= a[idxi][k]/pivot
      a[idxi][k]:= aop
      for j=k+1 to n
         a[idxi][j]:= a[idxi][j] - aop*a[idxk][j]
      next j
   next i
next k

for i=1 to n               "RHS transformation
   idxi:= 'idx[i]
   sum:= b[idxi]
   for j=1 to i-1
      sum:= sum - a[idxi][j]*x[j]
   next j
   x[i]:= sum
next i

for i=n to 1 step -1       " backward substitution
   idxi:= 'idx[i]
   sum:= x[i]
   for j=i+1 to n
      sum:= sum - a[idxi][j]*x[j]
   next j
   x[i]:= sum/a[idxi][i]
next i

call Mdelete(a)            " deleting auxiliary lists
call delete(b)
call delete('idx)
call delete('scale)
/E
E
#4
:n
gelpiv2
:z
" Notation: see the Procedure Function 'gelpiv1'

" Description: This procedure improves iteratively on the first solution
" from 'gelpiv1'.  The iterative cycles repeat until the maximum relative
" difference between two consecutive solutions is better than 1e-14.

call gelpiv1(n,a0,b0,a,b,x)
improve:
for i=1 to n
   'bb[i]:= b0[i] - dot(a0[i],x)
next i
call gelpiv1(n,a0,'bb,a,b,'xx)
flag:= 0
for i=1 to n
   flag:= max(flag,abs('xx[i]/x[i])>1e-14)
   x[i]:= x[i] + 'xx[i]
next i
if flag goto improve
call delete('bb)
call delete('xx)
/E
E
=r
#1
:r
"  ***  Linear  Equation  Solver  with  Iterative  Improvement  ***
:s
C
#3
:r
call gelpiv2(10,'a0,'b0,'a,'b,'x)
#5
:r
" This model solves a system of 10 linear equations for 10 unknowns by
:s
C
#6
:r
" Gaussian elimination with scaled partial pivoting and subsequent iterative
:s
C
#7
:r
" improvement.  The correct solution is  1,2,3,...,10; you may observe the
:s
C
#8
:r
" improvement of the solution (over that from simple call of 'gelpiv1')
:s
C
#9
:r
" in the bottom window.
:s
C
=p
=t
#1
:n
matrix
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
""
:u
#1
:l
a01
:w
6
:h
""
#2
:l
a02
:w
6
:h
""
#3
:l
a03
:w
6
:h
""
#4
:l
a04
:w
6
:h
""
#5
:l
a05
:w
6
:h
""
#6
:l
a06
:w
6
:h
""
#7
:l
a07
:w
6
:h
""
#8
:l
a08
:w
6
:h
""
#9
:l
a09
:w
6
:h
""
#10
:l
a010
:w
6
:h
""
E
#2
:n
auxmatrix
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
""
:u
#1
:l
a1
:w
6
:h
""
#2
:l
a2
:w
6
:h
""
#3
:l
a3
:w
6
:h
""
#4
:l
a4
:w
6
:h
""
#5
:l
a5
:w
6
:h
""
#6
:l
a6
:w
6
:h
""
#7
:l
a7
:w
6
:h
""
#8
:l
a8
:w
6
:h
""
#9
:l
a9
:w
6
:h
""
#10
:l
a10
:w
6
:h
""
#11
:w
6
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
%Ol,f,13,0,0,1,9>f,0,1,4,1,0
